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DINESH K. KAUSHIK 



ABSTRACT 


Ths burning rate o-P a comacsiie solid oropellant is affected be tbe 
cross— floi'. of oombusticn gases. This Dnanori-.encT, is called erosive burniAc,. ir 
Ihis .iork, erosive burning of composite solid pr-ooellant is in.-estigated be a 
steady, two - dimensional, chemically reacting, lurbuler-t boundarc, layer c -er a 
propellant surface The tu'^O'-'Ience is Tioaelled through the much used ti-io 
equation ‘'k - model while bpaiainy s Eddy Break - up model is employed to 
model tne gas phase reaction Jpodel . The method used to solve the boundary layer 
eauations is due to Patankar and Spalding, This is a semi - implicit -finite 
difference scheme in which the corservation equations '‘trans'f'ormed in x ~ u 
coordinates) are integrated over the control volumes constilutirg the boundary 
layer while marching in x - direction. 

The increase in free stream velocity brings the location of peak heat 
release zone closer to the propellant surface, thus, increasing its burning rale 
Also, the effect of normal burning rate, particle size and ,surf3ce roughness 
height is investigated. 
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= Mass fraction of the fuel in the composite 
solid propellant 

= Maas fraction of the oxidizer in the composite 
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= Stream function 
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Introduction 


Chapter 1 


(l.a) General 

When 5 solid prepellant burns in the presence of crossfloi’i of 
combustions gases, its burning rate is affected In the absence of a crossflow, 
bur'-'ing rate of a composite propellant depends on the chamber pressure, initial 
propellant temperature, fuel to oxidiser ratio, propellant configuration and 
oxidiser particle size. But in the presence of a crossflow, after a tl'ireshold 
velocity is reached, burning rate becomes a strong function of '-'elocity ?; 
usually increases with velocity (other paramers remaining fixed) This 
phenomenon is called erosive burning. 

High performance rockets require high thrust and short burning 
time which in twin require high volumetric loading fractions High loading 
fractions result in low part to length ratios and this means that crossflow 
velocity becomes high enough resulting in erosive burning. This may give rise to 
unequal propellant-web burnout, early exposure of the part of the rocket motor 
casing to the hot combustion products and overpressurization immediately 
following ignition. So erosive burning effect needs to be accounted seriously 
shile designing the grain and rocket motor. The phenomenon of erosive burning 
has also received a considerable attention because of nozzless rocket motors 
(which offer significant economic advantage & simplicity over conventional 
motors). Here gases reach sonic & supersonic velocities over the grain surface 
leading to high erosive burning rates. 


i 




(l.b) Background of the Problem 


As repored by Von Braun, erosive burning studies of solid 
propeiiants were made in rnany countries like Germany, Jspan, Brit sin, USSR etc 
and date back as early as iS39 The erosive burning effect was observed by 
Marsel on cordite. Muraour rias probably the first investigator to describe the 
erosive phenomenon explicity (Ref. (i)). 

Uver the years, a large number of erosive burnirrg models ha‘/e 
been reported. Generally these models differ itn the physical basis attributed as 
a reason for enhancement of the burning rate with increase in '.'elocitu 
Exhaustive literature reviews on erosive burning problem are published by King- 
and ‘<uo and Razdan^ Some of the recent models relevant tor the present work 
are reviewed I'lere in brief. 

Lenoir and Robillard^ developed their model on the premise that 
erosive burning is driven by increased heat transfer fromi the mainstream gas 
flow associated with increased* heat transfer coefficient with increased mass 
flux parallel to the grain surface. They calculated the total burning rate (r^,) as 
the sum of the normal burning rate (r^o^ and a second erosive term resulting 
from heat transfer from the core flow to the propellant surface i e 
Tb = rbo+ Tfee 

In this aquation it is presumed that the pressure-dependent normal burning rate 
rj^jj is unaffected by an increase in total rate at a given pressure. Also t-fce is 
postulated to be proportional to heat transfer coefficient This model has been 
widely used due to its simiplicity but the additive assumption 'll) is disputed 



and 15 5till not resolved. 


A model of composite propellant combustion was pi esenled by 
Hermance'’ He considers the combustion problem on a sutficientlu broad scale 
to be free of empirically derived constants. But he assumed a very unrealistic 
assumption about propellant surface a heterogeneous reaction is assumed to 
occur in a cusp between the oxidiser crystals and the binder, the cusp resulting 
from the spherical regression of the oxidiser particles Recent studies of 
extinguished surfaces by Boggs and others have shown these assumptions to be 
invalid 

Another somewhat more advanced combustion modal was put forward 
by Beckstead, Derr and Price^ This model is .based on a flame structure 
surrounding individual oxidizer crystals Three separate flame zones are 
considered 

(1) .A primary flame between the decomposition products of tthe binder and tbie 
oxidizer 

(2) A premixed oxidiser flame 

fS) A final diffusion flame between the products of the other two flames 
The oxidizer decomposition is taken as being the over-all controllilng factor in 
the combustion process. The results obtained matched e;^cellentally with 
experimental trends in the case of surface temperature and the effect of 
oxidizer concentration. But the predicted effect of particle size is somewhat 
greater than observed experimentally 


King^ has based his model of erosive burning on the convective 
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bending of the diffusion-controlled flsma zone As the total combustion length 
scale IS assumed constant, the resulting flame zone lies closer to the ui opellant 
surface and provides increased heat transfer He has criticized the Lenoir- 
Robiliard theory in the light of experimental evidence, as L-P model predicts a 
strong dependence of burning rata on the cere gas terriperature But King's flame 
bending model is found to underpredict tha-erosive burning effects severely 

In Lengelle's modal®, the basic propellant combustion mechanism is 
the granular diffusion model in which pockets of fuel vapour lea^'e the surface 
and burn away in an oxidizer continuum at a rate strongly dependent upon the 
rate of micro mixing of the oxidiser vapour into the fuel vapour pocket The 
effect of crossflow is assumed to be associated with the increase in turbulence 
with increasing crossflow raising the turbulent diffusivity in the mining region 
(thus increasing rata of mixing of binder and oxidiser product gases) and 
raising the effective turbulent thermal conductivity, which, in turn, increases 
the heat transfer rate from the flame to the surface Lengelle has ignored the 
ammonium perchlorate monopropellant flame. Also in his treatment of boundary 
layer, he has used the power velocity law all the way from the treestream to 
the surface on the groound that since this law is valid upto "inner zorie" of the 
boundary layer and gradients are steeper near the surface, this velocity law 
can be extended right upto the freetream. But this is an oversimplification of 
the problem in a bid to get a simple final expression for the burning rate. 

f 

The objective of Beddini's model^ is to investigate the role of 


turbulence in erosive burning. His flow field analysis is coupled with a simple 
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model of propellant combustions in which the mass burning rate is assunied to 
be di'-'ectly proportional to the heat flux from the gas to the surface and thie 
gas phase reaction is assumed to be described as a single-step homogeneous 
reaction, with no consideration of the heterogeity of flame structure associated 
with composite propellants 

Kuo and Razdan® used a second order closure model for 
characterisation of the flow field in erosive burning situations , the closure 
model ( k-e, model ) being used differing from that being used by Eeddini They 
have used Spalding’s Eddy Breakup Model for rriodeling the gas phase chemical 
reactions They concluded that propellants with lower normal burning rates are 
more sensitive to erosive burining and rate of burning increases with surface 
roughness 

Renie and Osborn^ presented a mathematical model for determining 
the erosive burning rate. The physical reason attributed to erosive burning 
effect IS same as of Lengelle that this phenomenon is due to an increase in the 
value of the transport properties in the gas phase reaction zone caused by the 
presence of a turbulent boundary layer. They have coupled a unique statistical 
combustion model, the petite ensemble model (PEM). They have also incorporated a 
turbulence model based on mixing length hypothesis as modified for surface 
roughness by Cebeci & Chang^'* . They have predicted the effects of cross-flow 
velocity, base burning rate, pressure and surface roughness on the erosive 
burning sensitivity. 

One recent model is due to Godon, Duterque and Lengelle They 



ha'.'5 aolved h-lavier-Stokes equations numerically, simplified with the shear flow 
appro'. irnation and for a two dimensional flow Mean the wall, they have applied 
couette flow approximation. They also use a mixing length closure model. In the 
"wall zone", they have obtained a flame height criterion by considering that the 
combustion is limited essentially by the diffusion between oxidising and fuel 
gases Both the flow description and erosive burning are sucessfully compared 
with experimental results cold flow-porous wall simulation with laser anemometry 
measured '-elocity profiles and propellants erosive burning measured with an 
ultrasonic transducer on a separated generator sample device They have very 
clearly predicted the effect of particle size, existence of threshold velocity 
and effect of normal burning rate. 

A very recent model presented by Sabnis, Gibeling and McDonald^^ 
employs a multi dimensional implicit Navier-Stokes. analysis. Most of the models 
which we have discussed used boundary-layer approximation to describe the flow 
field In a number of instances, such as the flow in the region of inhibiters. 
slots and submerged nozzles, the boundary layer approximations are not valid In 
the present model, ensemble-averaged Navier-Stokes equations for a cylindrical 
coordinate system are applied to the simulation of the steady mean floi-) in 
rocket motor chambers. They have used a two equation turbulence ( k-e ) model. 
The computed results show good agreement with the experimental data from cold- 
flow experiments on the axial velocity profiles Also the axial location of 
transition to turbulent flow is predicted reasonably well. 



(i.c) Physical Basis of Erosive BLrrdng Pheromenon 


While reviewing the work of different researcher in the previous 
section, we noticed three different reasons responsible for enhasncementt of 
burning rate with crooss-stream velocity 

(i) Heat transfer from a "core" gas to the propellant surface by forced 
con'.'ection 

<■ 11 ' Convective bending of the diffusion-controlled flame zone. 

(Ill) Enhancement of transport properties in the gas phase reaction zone by 
cross-flot'i induced turbulence. 

The first factor , if taken as a basis for erosive burning 
phenomenon, predicts a strong depends of erosive burning rate on the core gas 
temperature which is widely criticized in literature (Ref 2,7 & 9). 

The second factor underpredicts the burning rata severely (see 
Ref .2,' , therefore, can't be accepted as a reasonable basis for erosive burning 
problem 

The third factor, in fact, is the most agreed upon basis for erosive 
burning phenomenon and many researchers (see Ref .6-1 i) have developed their 
models with this as a reason for erosive effects. Kuo and Razdan® quite clearly 
showed that turbulence does affect the burning rate quite strongly Increase in 
free-stream velocity brings the location of the peak turbulent intensity closer 
to the propellant surface, which means that the turbulence eddies with high 
frequencies come closer to the propellant surface. This causes an increase in 


the mixing rate of oxidiser and fuel species, therefore, increasing the gas phase 



reaction rate 


In thie work also, this is the third factor, which is considred the 
ohusical basis for erosive burning phenomenon. 

(i d) Modelilng Approach 

As already stated, erosive burning phenomenon is observed as a 
result of the interaction between fluid freestream velocity and the combustion 
process. Any attempt towards erosive burning modeling should be started with 
integrating a fluid mechanical analysis with a combustion mechanism. In fact, it 
IS the surface temperature which provides the link between the burning rate and 
the gas dynamics In this work also, the approach is similar and we can 
articulate the calculation procedure in the following words- (i) calculate heat 
flux from gas to solid phase by showing the gas-phase conservation equations Sj 
thus determine the propellant surface temperature and (ii) utilize the knowledge 
of su'^'face temperature to calculate the burning rate. 

An acceptable model must account for the physically observed facts 
about erosive burning phenomenon in <i) a trendwise way (iil quantitative way. A 
few common observations of importance from the past experimental studies are 
enlisted below 

(i) Plots of burning rate versus gas velocity or mass flux at constant pressure 
are not fitted best by a straight line. 

net 

(li)There is a threshold velocity below which erosive effects are^observed. Also 
higher the pressure, lower is the value of threshold velocity. 
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(ill/ Siuwar burning propellants are more strongly affected by crossflows than 
higher burning propellants 

(iv' Small particle size AP monopropellants are less affected by the cross flow 
while the large particle size are affected tremendously. 

(v) At high velocities, erosive effects enhance with increase in surface 
roughness 

In this work, objective has been (1) to calculate erosive burning 
rate from an aerothermochemical analysis for a composite propellant (2) to 
examine whether the theoretical model meets the test of observation or not. by 
comparing the results with the above-stated facts from experimental studies. 



Chapter 2 Formulation Of The problem 


(2 a) Description of the Model 

In this work, an aerothermochenracal analysis to erosive burning 
problero IS presented in terms of two equation model t k - e > for turbulence 
and scalding's eddy break-up (EBU)^® model for combustion process, which 
considers the heat, mass and momentum transfer in a chemically reacting 
turbulent boundary layer. 

Hot gases flow over two dimensional slab of propellant and form a 
boundary layer on it (Fig. 1). To analyse the flowfield, where the heat and mass 
transfer is taking place between the boundary layer and the propellant surface, 
we make the following assumptions' 

fi) The boundary layer is two dimensional, quasisteady and chemically reacting 
(ii) The flowfield is statistically stationary i a. each instantaneous variable can 
be replaced by its mean and fluctuating part. 

(ill) Body forces are absent 

''iv) There is no radiative heat transfer. 

(v? The Lewis number is unity. 

'"vi^ *^ick’s law of diffusion is valid. 
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(vii) 


dU 

3y 


a ^ 
^ 3y 


« i 


andl^ 
jO u'v' 


« i 


In each case, these ratios can be assumed to be less than 5% for Mach numbers 
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less th5n 5. if the correlations between viscosity and velocity gradient and 
bett'jeen density and u’v' fluctuations are at the most 0 3, a conservative 
estimate ‘Fef 19, p 73). 

<viii'i There is no reaction-generated turbulence. 


Starting with the general conservation equations for a reacting 
compressible flow, we get the following simplified form of conservation 
equations, by doing an order of magnitude analysis combined with the above set 
of assumptions (See Ref 18, Ch.9 and Ref. 19, Ch 3 for details). 


(i) Continuity Equation 


~ ( 7) u > + ~ p ^ "i =0 

dA' oy 

(ii) X-Momentum Equation : 


u + p V 


(iii) Species conservation Equation -• 

’eff ciy 

where , i .= 9, for fuel 


- _ dY - .. dY 
P u -f- P V ^ 

oX dy 


ay 


(fe) 


iy ] 

‘ ] + ii. 


■ 1 ) 


( 2 - 2 ) 


0 , for oxidiser 


(2-3a.) 

rZ-Sb) 


But the species consevation equations are solved forYp and W lEq.2.3b) where 

V _ V V 

Vf 

This choice of variable eliminates the non-linear term in the equation for this 


variable. Also no separate conservation equation equation is needed for the 
product species as Yp j=(i - Yp - Y<> ). 


(iv) Energy Equation ; 

T n 3H , 7i .r afi _ a 


^ ( FF L, §5 + [ - ( ^ ). 

i, ^ 


l4] 

Pr 9y f 


(2-4) 
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(v) Equation of State : 

- _ 5 R" T 
~ l‘J 


(2 >5) 


For '.‘arioua correlation terms, we have used the following correlations 

d; PeynolcJs Stress u'v' is related to turbulent viscosity by the following 

equation 

- 7j U7- = Ait P t2-6) 

3y 

Thus = /i + 


(ii) Mass transfer by turbulent diffusion , that is , if) \jy Y,‘ is modelled in the 
follot*jing way 

- (puY Y = ^ ^ i'2-7) 

Sct ay 


dll) Heat transfer , (fi v)' h- is modelled by using Reynold's analogy , 

Ait ^h/ 


ipi-'Y h' = 


Prt ay 


n't 

U*o} 


(2 b) Turbulence Model 

For closure of the above stated problem , a two equation model of 
turbulence is used in this work . This model characterises turbulence by 
turbulent kinetic . energy ( k = ) and its dissipation rate 

(a = At ) The ( k - a ) model was originally developed by Harlow and 

P 

Nakayama fRef il) and was subsequently improved by many researchers Here , 
k-a model due to Launder and Spalding^® is utilized The potential advantage of 
such a transport model is, that it is capable of characterising a turbulent flow 
adequately without specification of the length scale required by mixing-length 


models 
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The turbulent kinetic energy equation for a steady , two- 

dimensional boundary layer flow can be written as ( see Ref 8 ) 

't Sk “ 


.0 U ^ 4- .0 


^ Bk _ d 

>J 1/ r — — ;r — 

3y oy 


( S ) I + ( i ) - ^ - 


lZ-5) 


The form of £.-equation^® is as applicable to uniform property flow . 
However , we modify it by retaining the full expression for the production of k , 
that IS . the second term on the right hand side of k-equation The e-equation 
for a steady , two-dimensional j boundary layer flow is (see Ref. 8) 

+ { §5 ) k“ ^ t 



Hhile solving Eqns. (2.9) and (2.10) , the influence of nearby wall is 
considered by means of widely used wall-function method (see Ref 16 for details) 
The turbulent viscosity Mt is related to k and a by the following relation , 


ilt — 



(2 'll) 


(2.c> Gas - Phase Reaction Model 

In the boundary layer above 
gases may react in several steps , but in the 
global forward reaction is assumed ; 

f- VqO -+ if^pP 

For such a reaction , the following Arrhenius 
reaction rate for fuel species is used , 



a burning composite propellant , 
present work ^ a single - step 


( 2 - 12 ) 

type expression for the global 

(2-13) 


where 


(2-14) 
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In species conservation equation <2 3) , we need to know average 
value of reaction rate , Up c uJq can be eliminated by the choice of variable , 
OF, as discussed earlier ) In this work Spalding's Eddy Break-Up ( EBU) model is 
utilized to derive an expression for . In this model , it is proposed that the 
gases in a turbulent flame , at high Reynolds numbers , should be considered as 
lumps or eddies of unburned and fully burned gases. The rate of burning 
depends the rate at which fragments of unburned gases are broken into still 
smaller fragments by the action of turbulence. 

The EBU concept can be used to model the gas - phase reaction 
rate ■*^or erosive burning problem The fuel and oxidiser gases are unmixed as 
they come out of the propellant surface .The presence of high lateral gradients 
in the boundary layer creates the turbulent eddies. It is reasonable to assume 
that these eddies engulf fuel and oxidiser gases , thus giving rise to fuel and 
oxidiser eddies 


Applying the EBU concept further , it turns out that in a diffusion 
controlled reaction , the rate of consumption of fuel is proportional to the 
dissipation rate of fuel containing eddies , characterized by 


I' 5 ■ k ^ g turbulent time scale. So , we can write , 


£ 

k 


Up oc — /O 

If we further assume that in the conservation equation for Y'^' 


(2.i5’i 
, only the 


production and dissipation terms are dominant , we have , 


a- ^ k j 

0 £ t cfy ^ 


(2 >16) 



iS 


Thi= a^Bumption is particularly good in in the near - rtall region This is the 
region ot importance because most of the chemical reaction tahes place near 
wall 

(2.d) Solid Propellant Burning Rale Equation 

We express burning rate , j of solid propellant by Arrhenius law 
of surface pyrolysis 

rt = AS exp ( - ^ ) , (2-17) 


(2.e) Boundary Conditions 


We need to specify the boundary conditions to solve the 
conservation equations at the solid - gas interface and freestream edge of the 
boundary layer. 

(1) At the solid - gas interface / wall 


^i> u = 0 


(2-i8a) 


^ 11 ) 



(ill) T = Toj 


(2- 18b) 
(2-180 


(iv) The mass balance equation at the interface can be written as follows 


total mass of the species k vapourised from the solid phase = mass of the 
species k transported away from the interface by the normal velocity in gas 
phase + mass of the species k .transported from solid to gas by diffusion 
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Or, 

( 7 I' V. )u 
whe'"s . 1 


* 2 D If L 


0 


(2-19) 


F or 0 


fv’' The energy balance is written as 


3Tp 

'W 


ay 


+ Ps '“fc ( hw- - h„+ ) 


(2 -20) 


In Eq (2 22) , first term, on the left hand side represents the net heat flux to the 
solid propellant . the first term on the right hand side represents the heat flux 
from gas to solid surface and the second term represents the net heat release 
at the surface. Subscripts ’ - ’ and ' + ' are respectively for the solid and gas 
side of the interface 

If we assume that heat conduction into the solid propellant is 
dominant in a direction normal to the burning surface , the heat conduction 


equation in the solid phase can be written as , 

^ = Ps Oj rt, for - oo < y ^ 0 

Sy*" ^ 

Integrating Eqi2.23) from y = -os , where Tn = T„j and =0 to y = □ , 

^ ~ ay 


( 2 - 21 ) 


aT, 


where Tp = Tpj , we get 

dip 




3y 


— /j; Cf rj, ( Tpj Tp^- ) 


( 2 - 22 ) 


The net heat release at the surface per unit mass which is equal to the 
diffe'-’ence between the enthalpies of the gas and solid fon their respective 
sides of the interface ) , can be expressed as ( see Ref 8 ) : 

Qs t Tpj ) = h«+ - hw- = Qj -f-' ( Cp - Cj ) ( Tpj - Tpj ) (2-23) 

where ,~Q^ is the net surface heat release at a surface temperature , Tpj . 
Combining Eqs (2.22), (2 24) and (2.25) , we get 


io 


■r t 




) T,5 ] 


(Z-24) 


(v'i) '‘he boundary conditions for k and e, are to be apolied near the wall and not 
at the wall because in the low Reynolds number zone in the vicinity of the wall . 
k - t 'fiodel IS not valid ( Ref .iS , p 684). 

Near a wall , the production and dissipation terms in the 

k -equation are equated , assuming that only these terms dominant ( a reasonable 

assumption at the edge of low Reynolds number turbulent region ) . Thus , 

£ = :^ ( |y f (2-25) 

p V 3y ) 

The eddy ‘Viscosity , near the wall is calculated using a modification of Von 
Driest formula , 


= ,3 [ Di, ( y + Ay ) ] 
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9y 


{2-26) 


The damping coefficient Di- , is altered to include the effect of surface 


roughness and the effect of local shear stress 

( y + Ay ) ]5 _T_ 1 

Tm ) 


Di, = 1 - exp ( - 


where , 


Ay = 09 


4^ - rJ exp ( - ^ ) j 


and , 


+ P Uj Rf; 


— 




{2-27) 


(2-28) 


(2-29) 


Eliminating from Eqs. (2.27) (2.28) and (2.11) , we get 


k = 


[ K Pi, ( y -f Ay ) 3 / ^ 


aI^ 


ray) 

I ay / 


(2-30) 



IS 


a = [ D., t y + dy ) ]^ 

(2-31) 

fll) At the freestream 


a'' u = Uc: 

(2 *323) 

(ii) T = T» 

^2-32b) 

Uii; = '<0 =0 

(,2 -320) 

fiv; = 0 

dy dU 

(2-32d) 


liith the above set of boundary conditions at the wall and 


freestream , the mathematical formulation of the problem is complete. Next we 
have to solve the set of conservation equations subjected to the above-stated 


boundary conditions 


I Ch apter 


Patankar - Spafding Transformation 


The set of conservation equations f 2.2 to 2 4 and 2 3 to 2 10 ) is a 
systeTi of partial difference equations which is to be solved sirnultaneously.This 
system is parabolic in nature and involves non-linear terms In this work , this 
system of equations is solved numerically by a very popular technique proposed 
by Fatankar and Spalding In this chapter, the conservation equations will be 
transformed to Patankar and Spalding cordinates. ' 

(3. a) The Choice of Co-ordinate System 

Fig .2 shows the boundary layer over a flat plate , and attached is a 
system of cartesian ( x ~ y ) coordinates The drawback of this coordinate 
system is that it puts only few points in within the boundary layer where it is 
thin and so can't provide high accuracy Ihere.Similar is the case with Von Mises 
coordinates ( x ~ ) , shown in Fig.S.What is wanted is , if''accuracy and 

computational economy are to be combined, is a grid which conforms better to 
the actual shape of the bopundary layer. For this , Patankar and Spalding have 
suggested a coordinate system in which cross-stream coordinate is represented 
by non-dimensional stream function , o , defined as 

u =r ' k (S-i) 

— 0 , 

where and 0 ^ are functions of x alone .Fig .4 shows this system ( x , tj ) of 
coordinates called Patankar and Spalding co-ordinate system To transform the 
conservation equations in ( x ~ o) ) coordinates , it is needed to transform them 
first in < X ~ i/f ) coordinates which is done in the next section. 
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BOUMDARt layer over a flat fl«te liith 
C00RDI!-!ATE3 



BOUNDARY LAYER OVER A FLAT PLATE WITH x 
CODRDII-4ATE3 




FIG. 4 BOUNDARY LAYER OVER A FLAT PLATE l-ilTH 

COORDINATES. 





(3 bi Transformation to Von Kises Co-ordinates 


The stream function Hi is defined by 

^ 0 '' 3 > 25 ) 

&y 




dX 


3 V 


with 


alli\ 

3v 


- li-ii 


3^ 


dx 


= - me 


<3 -lb) 


(3-33) 


(3-3ta’l 


We have to transform ,say a function f in ( x , y ) coordinates to ( x , u > 
coordinates.From Calculus , 



Dr (Simply in operator form 



Since f ^ ] “1 «nd making use of Eq (3.2b> , we get 

\ 3x 'y 



Similarly , making use of Eq.(3,2a) , we get 


With the help of transformation relatiorte (3.5) and (3.6) , we can 
.rile the ronsa^vation equations (2.2 to 2.4 and 2 9 to 2.10! in the follomns 
general form (droooing the subscripts in Eos. (3.5) and S 6) , asaun,ing that their 
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meaning ie understood ) 

^*4^ _ a j' 7 ) - p \ , i o 

I i + td 

where #' stands for any of the flow variables , is the source term 
corresponding to the variable and effective value of the transport 

property in the flux law of <p. In Table 1 , and for all the 

conser'.'ation equations are shown. 


Table 1 : Conservation Equations in ( x ~ ^ ) Co-ordinates 








(o.c) I ransi formation to Patarrfcap - Spalding Co-opdinate System 


Here . the ppoblem is to tpansform a function f (x , m to f <: 
I'la c=n 1 ‘jrite equations ( in the openatop foptn ) similar to Eq.G 4) 

^ dX Agi) ^ dX dX . dU) 

j^l =f^)('32£] 

^ d^/.' ^ cix 3u} \ 3aj /x^ 30 ix 

Since . ( ) = i and ( |^ ] = 0 

^ dx / ^ 3c*) f ^ 

the above equations reduce to 

i 3.] =( 1 .^ + 

> Bx /«} ^ 3x ■'(,} % 3ijj )x^ dx •z® 

i ) = (I; K i L 


To e’.'aluate 


dOJ 

and 

3u 

- 3x^ 


d0_ 


, we write Eq.O.l) as 


= ti»| + u) ( ) 

remembering that 0, and 0g are functions of x alone . 

From calculus , total change in as a function of x and u can be written as 

=■» = { i! I"* + ( E 1='“ 

If ij} IS constant , then d0 = .0. Therefore , from Eq.O.li) 


s , u). 


(3-8) 


(3-9> 


C3-i0) 


(3'ii) 



2S 


FrC'iM Lq ^ z ‘ ID' 


dK 


CvVi H 


and. 


( i, I = < ) 


Using Eqs *3 13) and (3 14) ^ Eq (3 12) can be expressed as 

il^) = ' 

• CiA 


( -IM 1 S + “1^ ^ W’E - 5^, ) 


( 3 - 13 ) 


G-i4) 


(3-15^ 


Further, from Eq= <3 ii) and (3 14) , we obtain 

I aju '| _ 1 _ 1 

I }y ( ^ 1 k — ^\ ) 

^ 3to Jx 

Makins uae of these results in Eqs. (3 8) and (3 9) , we get 

< I: ^ ^ ' (1^ I 

( 5. \ ^ 1 / 8_ ] 

' ot) .’x ( i/'E - i('i ) ' 'x 

where, 

. _ i d^/ii 


(3-16) 


{ 0E - 01 ) dx 


(3-17) 

i3-i‘8) 

i'3-19) 


( 0E - '/’I ) 


(3-20) 


( 0E — ^ 

With the help of the transformation relations (3 16) and (3 17) , the 
general form of the conservation equations , Eq.(3.7) , can be written in the 


following form 

n _t 

+ ( a + b u ) 


dX 


dCi) do 


a / d u Fji.eff 3^ 

I ^ — 1 


Or, 


Q4 j / = j_ K , 1 ^ 

- +(a +bca)5j^ 


A 

3u) 


( 0E — ) 




M ^ 4- J- 
3u iO 0 


(3-21) 



With 


( We - W, f 


( 


2 ) 


d 



(3-23) 


In T able 2 , c and d for all the conservation equations are presented 

Eq,(3.20) is the general form of conservation equations. Some 
remarkable features of this equation are ; 

''i) In view ot Eqs (3.3a) and <3 3b) , we have from Eqs. (3.iS) and (3 19) 


3l 


_ 

~ ( We - Wi ) 


(3-24) 


b 


nriE 


rn. 


i) 


( We - Wi) 

where the mass fluxes rh^ are tri, are the mass flow rates per unit area divided 
by ( We - W\ ) across the grid boundaries, and are positive for flow in positive y- 
direction.Thus , t a + b w ) denotes the mass flow rate across the line of 


constant w , divided by ((/»£- l(), ). 

Cii' Patankar - Spalding coordinate system has an extra advantage over the 
choice of (y/yg) as a variable (if chosen , instead of o). For a constant x , 
( a + b w ) is a function of u alone , but with (y/yg) as avariable , this 
coefficient would have involved an integral over that variable of the forward 
velocity , thereby , making this coefficient more troublesome to evaluate. 










Chapter 4 


Method Of Solution 


As we have seen in the previous chapter , all the conservation 
eciations can be represented by the general •form , Eq (3 20) .To solve this type of 
equation which is parabolic in nature , we use the methodology due to Patankar and 
Spalding This is a marching technique , that is , the finite di-fference eqations are 
sol'.. ed at the succesive x-stations, using the solutions at the previous x-location. 
In this chapter , first transformed conservation equations (Table 2) are cast into a 
set of finite difference equations and then some salient features of the solution 
procedure will be described in brief. 

(4 a) Finite Differencing of Equations 

In Patankar and Spalding method , implicit finite differencing is done 
which connects the value of a dependent variable , , at a downstream location 

node 1 , with the values at the two neighbouring nodes at the same x-location , 
■’^,- 1,0 ^' 1 + 1,0 3 ind with the values at the three corresponding upstream nodes , 

^ 1 + 1 , u- Equivalently, it can be written in the following form . 

where A, , B, , C; and D, are coefficients for which expressions are to be derived 
by integrating Eq.<3.20) over an appropriate control volume (see Fig.6)lt should be 
noted that the coefficient q carries the effect of the upstream (i.e. previous x- 
location) node's). 

Before we start integration over the control volume , it is 
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nece==3ry to employ a ^-profile assumption We show this assumed profile in Fig 5 
Hslfwiy between idjBcent pairs of grid points in the range 2 to fN-i) are drawn 
‘.'erticai broken lines , dividing the whole u)-range into (N-i) intervals The value of 
any '.'ariable ^ is supposed to be uniform within the interval. 


The control volumes , over which we are going to integrate 
Eq IGkare illustrated in Fig 6 The boundaries of all the control volumes except 


the two near tl'te edges are represented by 


upper edge . 

i 

0) 1 i ( W,, -f* 0). ) 

X. 

f4 2) 

lower edge • 

^ 1 = 1 ( U, + ) 

■4 3) 


The lower edges of the control volume near I boundary (i.e for i equal to 2) is 

w i = 0 ^4 4) 

‘'a 

Also. the upper edge of the control volume near the E boundary (i e for i equal to N 


-i) IS 


u) 1 = 1. 


(4.5) 


We •.-.•ill satisfy conservation equations for each of the control volumes and thus 
conservation will be satisfied for the whole boundary layer .Integrating Eq C3.20) 
over 1 *^^ control volume , we get 


i+l 


JL 

5.x 




term (1) 


+ (( a + b(o ) f 
^ v-iM 


term <ii) 



I boundary 


E boundar\ 


;.5 # - PROFILE ASSUMPTION FOR FINITE DIFFERENCING 


N 





THE CONTROL VOLUMES USED FOR THE DERI'.'ATIDIT OF THE 
FINITE DIFFERENCE EQUATIONS 






term (v) 


(4 & 


where subscripts U . D and M stands for upstream , downstream and 
intermediate :c-location.With the ^-profile assumption of Fig 5 , various terms of 
Eq,(4,6;' can be written as follows 



the .c-direction convection of ■ 



DD 


(i) + (iii) 




L ' 


E 


llLo . . f ‘ “^-1 ] 


In visi-i of Eq'3 23) and Eq (3.24) 


'■ a + b UJ ) = :tt — =^rT- [ tfi, + 03 ( rrig. - m, ) ] 




HO tl"‘at 


a + b o.'| ^ = 


"'.+i 


i+l (t^g; “ 


whei'H , fu 1 

1+7 


rii, + u) 1 ( mg: “ ^ 

'+2 


m. 1 


Similarly , a + b oil = 




With , m. 1 


ni, + o3_^ ( rrig - ffi, ) 


Therefore second term can be written as 


(ii) 


[ 




tt!ig - {i»,)y i+5 ' i+l 


1 I ' 
‘*2 ‘“*5 

Jm 


Choosing t’1 as the downstream location , that is , 


and , +■ ) 

4, i ^ ( ^,-i *+ ^; ) 

*"5 

With these i ^4 15) can be written as 

(ii) : ^y:. ^ - - r - f m, . 1 + ( tit 1 - m i ) #■; - m i #i-i ~| 

- ^Vu j _ '+2 '+2 '2 ‘2 Jm 


(4 .S'* 


(4 10 3 


^4 li> 


(.4.12) 


f4 13) 


(4 14) 


(4 15) 


(4 16 a) 
(4.16b) 


(4 17) 
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j(«,0 ciirfuHi'ye term (vA can be rewritten as follows 


" r i£ " 


i 

r 

iwd/ 

dW 

l+|,M 


^ 3y 


?nn 


&to 




- 0,?i 


lU 


^,eff 3y 


|+|,M 




(4 IBl 


<■4.135 


If it 1= asHurned that . 


'Ci0 

&y 


d£ 

dy 






< >D 

< yi+i - y. 'u 

< #; - ^i-1 

( 'Ji - ‘Ji-l 


(4 20a) 


14 20b5 


term U'/') csn be written as 


(M : 


±. 




T,^ ( >D - - ^i-1 '^D ] 


fA 21) 


U| 

■^.+1 ' < y^+x - y; t 


i4 22) 


‘ I^,eff li.t ■ t4,23) 

' "’’i-l ' ( y. - ys-i t 

The source term (v) . represents the effect of sources of # within the 

t. r^£ V Althnuah sources are non-linear functions 
control volume par unit increment o 


of in general , these terms are 


evaluated by adding in upstream values a term 



dspc'^'diri'S lititiHsrly on ip 


Ci;,o d(j 


= Si + s; ^.,c 


(4 24^ 


In the next section , source linearization procedure i-Jill be presented 
, that 1= . the expressions for S; and S/ for each individual conservation equation 
will be derived 

Moh , with the help of (4.9) , (4.17) ^ (4 21) and (4 23) . Eq.(4 6) can be 
written in the following form- 


, ( u) , - u» , ) 

( - V'l^O ’"i , i , ^ s 

1 'II'e - b,),j £x Z'd'E - 0|),J <+2 >-2 


+7n — '' » ' V" ^ T., 1 - T i ) - 5 ' > y Q 
- (^|)y 1+2 l"2 j 


( rr — S tt* ‘ T 1 - m x ) \ ^i+i,D ^ ^ ^i-i ^ 

t ttilg. - 1+2 1+2 > ' t We - Wo ' z ' 2 '' 


f W , - U) I ) 

‘+2 * , c 

K ^i,U * ®i 


(4 25) 


The first two terms in the coefficient of can be simplified with the help of 

Eqs <3. 18), (3. 19), (3 .23) and (3.24) in the following way 


we have , 


We * Wo 

((i^E “ 


1 . ^ --!!! sx 

^ (W - ^’l>U 
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< W i - U) jL ) 


ir.£r£for£ , — Li L2__ 

"■^'e - ''Vu 


6 x 


( w i - u) , ) 


-( m ,1 - rh 


(4 26) 


Sx ZC^E- - ii|)y "'.+1 ^ 

llith U-e h5l?j of Eq '4 26) , Eq (4 25) can be written as (after multiplying both sides 


by C.-. ) 


’2 


< T^i - rh .1 ) + 


8x 


.+ % ^(('e ((‘Vu '■^2 '"'+1 ' ‘ - 


( T_i + rh 
‘■*2 


=-=4 




-{« 


1 T 1 - rh 


Vi>} <^(+1,0+ ( 


Sx 


‘+1 ’‘+1 ' ^ ' ’ U0E - ^^,>u 


< T. 1 -f m. 1 ) } 


I-l/D 


I ' “.i ^ + s, fix I 


(4.27' 


Comparing Eq<4 27) with Eq(4.i), expressions for the coefficients A, , B, , C, and D, 
are obtained 


«■ 


§2 ( T, 1 

li^g. - ’+2 



(4 28) 


B 


— — ( T 1 
((('e - i -5 


m, 1 ) 

1-5 


(4 29) 


C, - ( U) , - u) i ) + Sj Sx 

i+g * 

D, = (u)j_-uJi) + Aj + Bj- Sj'Sx *‘4 

i+- i-g 

The coefficients A. .and B; are to be modified as these may become 
negative for high lateral mass flow rates and thus , may yield physically 
implausible values of ^^ 1,0 .The following modification is incorporated in this work ^ 
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Al=o Irsnsfer coefficients T 's are to be modified at the boundaries 

becsuse 

■ i'' 5 f E-L":''ir.d’!r'j , tt-.e snd adjoins the undisturbed flow and ^i^-gradients are zero 
(ii) st I-boundsrg , due to the presence of a wall,there are sharp variations in 
transport properties. Also, special effects such as kinetic heating or chemical 
reaction , make the diffusive flux proportional not to ( - 4>2 '* but to some other- 

difference 


So, for 1 = 2 


^-| ” i 


(4.34) 


and , for i = M-l ■ T 


.J.1 

'+2 


0 


(4 35) 


where J I .-tunds for the diffusive flux of 4> at the 1 boundary and represents 

fl, 

the above stated wall effacts.For the E-boundary , ^4^ is zero. 

(4 b) Source! Linearization 

The source term 01,0 can be expanded by Taylor senes i upto second 
terro only as it is assumed that source of * depends linearly on f ) in the 

following way 
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rhen we can obtain the following expressions for the coefficients S, and b. 
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In Tablfe 3 , c*, and S| are presented for various conservation equations 



where , F 
and , 
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For fuel species conservation 
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equation , source 


linearisation is done in the 


folloHing way 
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or 
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i + i , □ 
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and =ub=equently tl-.e coefficients A; and B, are modified for this equation. 


(4.C) Solution Procedure 

Tl'ie finite difference equation (4.27)' gives rise a system of 
algebraic equations which are tri-diagonal in nature. These are solved by the 
tri-diagonal matri.': algorithm (TDMA) described in Ref. 2i. Around eighty grid 
points i-jsre considered and half of them were distributed with in iQ% of the 
boundary lay«i' thickness, where the flow variables have sharp gradients. The 
forward step-size is calculated in two ways; (i) based on boundary layer 
thick'-iess (equal to 0 3 times the boundary layer thickness ) (ii) based on the 
mass floii rate at ! Ss E boundaries (Ref ,20) . In the program the actual step size 

was taken to be minimum of the two. The initial wall values of temperature, fuel 

\ 

and oxidation mass fractions were guessed and then iterations were performed 
until the boundary conditions at the wall are satisfied 

In Arretvfi<-C, the flow chart and descriptions of various 


subroutines is presented. 



Chapter 5 


RcsL!!ts,Dlscussions And Conclusions 


(5 a) Rtasults and Discussions 

Mutiierical computations were done for two types of composite 
propellant; compositions <D Ammonium perchlorate (75'/,) and PBAA/EFON (25%) and 
(2) A,nf..ora'ji.i perct'.lorate (65*/.) and Polysulphide (35’/.) The basic aim behind these 
calculations is to predict the effect of freestream velocity on the burning rate 
unde-- tyi.-'icul flow conditions, similar to those prevailing in an actual rocket 
motor This was done by varying the freestream velocity in a reasonable range 
The arious physical properties which were used in calculations alongwith tlie 
constants used in turbulence modeling are listed in Appendix B The details of 
the t rccedures followed in obtaining some of the properties are described in 
Pef 8 While calculating the paraineters associated with the glotail swingle -t-p 
fon., 5 ,'d r.iSCU0,Mt ll Msumsd lhal primary flame collapsee at the propellant 
eur(ac=i«rely dapoa.lina ita heat there ae the gasifioation takes placeisee 

R«f 

{ ton-dimensional velocity Cu/Uoo) and tampei ature (T Tpg- ca- ps' 
profiles are shown in Fig.7 and Fig.8. The velocity gradient is ma..im 
propellant surface and decreases monotonically to nedily o ^^alue a. 
freestream edge. The temperature rises rapidly in the near-wall reyiot. 
value at the propellant surface and then gradually approaches the freestream 
value Tine rapid rise in the near-wall region is brought dbout by th_ 
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re=ctioni tjl *rv_, dace timers Fig.9 shows the calculate! distributicna of fuel and 
oxidiaer r„.5a r-ractionalt is evident from Fig 9 that fuel and oxidiser mass 
fractions drop to '.'ery small values near the propellant surface implying that 
chemical reactions are quite fast in that region and therefore.the rate of heat 
generation di.ie to chiemical reactions in gas phase also rises sharply from its 
zero jslue at the surface , becoming maximum in the near-wall region and then 
gradually decreasing to zero at the freestream edge of the boundary layerfthis 
15 quite consistent with the temperature distribution shown in Fig,3). 

It is also clear from Fig 8 that temperature gradient becomes 
steeper as tl'ie freestream velocity increases which means higher gas-to-solid 
heat flu.*' This fact is expected because higher freestream velocity induces 
higher turbulence intensity near the propellant surfaca.This causes an increase 
in th'e iia..irig rata of the oxidiser and fuel species and therefore gas-phase 
reaction rate increasaa.Thus we have a double-fold effect of increase in 
e;-:ii .elocity firstly peak value of rate of heat generation due to 
chemical reaction incraasas and secondly the location of this peak value comes 
closer to the propellant surface.Both of these factors add to the gas-to-solid 
heat flux md th.us makes the temperature gradient steeper in the near-v-jall 
region. 

Another ooriseauence of increase in freestream velocity is the 
higher total burning rateCsee Fig.iO) or increase in the erosive burning 
augmentation factor <rj,/rj,o) as shown in Fig. 12 .As discussed above the gas-to- 
solid heat flux rises with the freestream velocity which means that propellant 
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=urr=ce feruoersture iticreases^making its burning rate higher On extrapolation, a 
threihOiJ elocitu ^belotM which there is no increase in burning rate) nearly 
eciuai f-O zoo nv's can also be noticed in Pig 12. 

Fig 11 shows the variation in burning rate along the propellant 
Sur's'-fi obaei ved that the burning rate decreases as one rAoves 

dowi-etresm This is brought about by the fact that the boundary layer is thin 
near the upstream end resulting in maximum gas-to-solid heat ^luxA 5 one moves 
doi/ji'istreaiu , boundary layer thickness increases and reducing the gas-to-solid 
heat flux 

in Fig 13 , the effect of surface roughness on the erosive burningi 
augmentation factor for different freestream velocities is shown As the 
roughness Iteight increases , the augmentation factor also increases This is due 
the fsct that increased roughness height causes more turbulence activity near 
thie propellant surface and therefore more mixing of fuel and oxidiser 
specie* Alsu the effect of roughness height is more pronounced at the higher 
freestream ,'elocitias.This is due the fact that at low freestream velocity , 
thickness of laminar viscous sublayer is more and submerges the roughness 
element But at highier freestream velocities , goes on becoming thin and 
roughness facor greatly influences the flowfield near the propellant surface. 

Fig 14 shows the effect of normal burning rate (r^^^ erosive 

burning augmentation factor. The curve with normal burning rate equal to 0 69 
offl/s is for AP/Polysulphide and the other curve with normal burning rate equal 
to 0.8 cm-'s IS for AP.'PBAA/EPON propellant.It is observed that the propellant 
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the - 


....we, r-,or.v.al burni ,-.3 rate is more sensitive to the erosive burning effect 
Mv:- p. c,'. ell ant with higher normal burning rate.This can be attributed to 
.->:t 3 fast burning propellant gives rise to thicker boundary layer over its 


su' '-.ue to liigl iSp '.rossstream mass flow rate than would be obtained in the 
case rr a slou burning propellant .In the thicker boundary layer, there is less 
^ei.'ipvr ature gradient and lower gas^to-solid heat flux which results in 
lows--" wroci /e burning rata. 


rtlso based on above discussion it can also be argued that as small 
parti-le Size AP composite propellants have high normal burning rates , 
therefore they are less affected by crossflow. 


(5 b) Cirwlu&ims 

It i£ observed from the present calculations that freestraam 
’.'elwcity greatly influences the burning rate of a composite propellant .Increase 
in velocity increases the gas phase . reaction rate and also brings 

the '«j\.»ion aone closer to the propellant surface .The effect of roughness 
h&ig‘*.t la to increase the burning rate and this effect is more pronounced at 
higher -‘j ■ h ••.-u'l velocities .The propellants with lower normal burning rates 
are wore sensitive to the erosive burning effect than the propellants with 
higher normal burning rates. 

The present work can be easily extended to take into account the 
effect pressure gradient .Also the work can be modified for axisymmetric 
bcvr diry l^yer . thit is . to model the erosive burning problem of cylindrical 


A I ■ '3 


c-ji.cc'ixie solid propellant .With the availability of more information about the 
cl f-t- .ic :•! hinetica of the gas phase reactions > there is ample scope for modeling 
tl.-s: * eaction rate in a better way. For example, a combustion model based on 
mult ipiu flsfiies can be employed to have more accurate estimate of burning rate 
ci ■ , j'j jite propellant (see Ref 5). 
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Appendix B 


Physical Constants 


[he .aiue= of the physical constant used are indicated below 


A5 

= 

26174 


Cp 

= 

03 

kcal / kg-K 

*■*^5 

= 

0 38 

kcal / kg-K 

Eas 


30 

kcal / mole 

P 

= 

iOO 

atmosphere 

Pr 


= 0.7078 

Prt 

= 

0.9 


Qs 

= 

-166 

kcal / kg 

Ru 

ss 

1.98 

kcal / kmole 

Sc 

s; 

0.7078 


Sct 

ss 

0.9 


Too 

« 

2250 

K 

Tp. 

St 

298 

K 

S 

* f&S 

ss 

800 

K 


M 

30 

kg / mole 

•Wo 

s 

27.8 

kg / mole 

t^fr 

ss 

1 


i£/0 

ss 

32266 




1600 

kg / m® 

Ci 

ss 

iD 


Ca 

sa 

1.3 
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i.57 


C3 

Ca = 2.0 


*--u) 


O.iS 


% = 0,03 

= i.o 



Appendix C 


Flow Diagram And Its Description 


15 =hows! the flow diagram of the solution procedure Hhat various 
Subroutines do, is described below in brief. 

MARCH.START 

I'i) Reads the input data 

Qi) Sets the initial profiles of flow variables except k and c 
PHV_PPOPERTIES . 

ti) Sets the initial profiles of k and e 

'll) Calculates the value of various physical properties, hi , , p etc. 

at all the grid points. 

EMTRAINMENT_RATE 

(i) Mass flow rates at I and E boundaries are calculated 
PEADV . 

<i> Calculates '/-distances between the grid points 
'll) Sets the step size, if next step is to be taken 
''iii) Sets the initial (guess) wall values at the downstream location 
WALL.REGION 

(i) Wall functions are calculated. 

SOURCE. 

(i) Calculates the value of Sj and 5) for all the conservation equations 
and at each grid point 

SOL*/E . 

(i) Tri-diagonal system of algebraic system of equation is solved. 
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8 CUKL-REGIDt -I 


(i) The coefficients A;, B„ etc, for each conservation equation are 
calculated 


^ii) Calls the WALL-REGION, SOURCE and SOLVE to give the value of flow 
variable at the downstream location 


3 CO! i'-'EPGENCE ' 


(i) The wall boundary condition are checked and iteration is done till 
they are satisfied. 
iO STEP-OUTPUT ; 

<i'> Unites the output of the program at a x-location 





FIG. 15 FLOW 


diagram of the solution procedure 
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